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Abstract. Identifying the processes that determine the initial mass function of stars (IMF) is a fundamental problem in star 
formation theory. One of the major uncertainties is the exact chemical state of the star forming gas and its influence on the 
] dynamical evolution. Most simulations of star forming clusters use an isothermal equation of state (EOS). However, this might 

be an oversimplification given the complex interplay between heating and cooling processes in molecular clouds. Theoretical 
predictions and observations suggest that the effective polytropic exponent y in the EOS varies with density. We address these 
, issues and study the effect of a piecewise polytropic EOS on the formation of stellar clusters in turbulent, self-gravitating 

molecular clouds using three-dimensional, smoothed particle hydrodynamics simulations. In these simulations stars form via a 
process we call gravoturbulent fragmentation, i.e., gravitational fragmentation of turbulent gas. To approximate the results of 
published predictions of the thermal behavior of collapsing clouds, we increase the polytropic exponent y from 0.7 to 1.1 at a 
critical density n c , which we estimated to be 2.5 x 10 5 cm 4 . The change of thermodynamic state at n c selects a characteristic 
mass scale for fragmentation M ch , which we relate to the peak of the observed IMF. A simple scaling argument based on 
the Jeans mass M s at the critical density « c leads to M ch oc n^ - 95 . We perform simulations with 4.3 x 10 4 cnr 3 < n c < 
CIh' 4.3 x 10 7 cirT 3 to test this scaling argument. Our simulations qualitatively support this hypothesis, but we find a weaker density 

dependence of M ch oc n -0- 5±0A _ We also investigate the influence of additional environmental parameters on the IMF. We consider 
variations in the turbulent driving scheme, and consistently find M s is decreasing with increasing ;i c . Our investigation generally 
supports the idea that the distribution of stellar masses depends mainly on the thermodynamic state of the star-forming gas. The 
thermodynamic state of interstellar gas is a result of the balance between heating and cooling processes, which in turn are 
determined by fundamental atomic and molecular physics and by chemical abundances. Given the abundances, the derivation 
of a characteristic stellar mass can thus be based on universal quantities and constants. 

5^ ' Key words, stars: formation - methods: numerical - hydrodynamics - turbulence - equation of state - ISM: clouds 
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1. Introduction 

One of the fundamental unsolved problems in astronomy is the 
origin of the stellar mass spectrum, the so-called initial mass 
function (IMF). Observations suggest that there is a charac- 
teristic mass for stars in the solar vicinity. The IMF peaks at 
this characteristic mass which is typically a few tenths of a 
solar mass. The IMF has a nearly power-law form for larger 
masses and declines rapidly towar ds smaller masses {Scalo 
1 1 998l lKroupal2002l IChabrie J2003I) . 

Although the IMF has been derived from vastly different re- 
gions, from the solar vicinity to dense clusters of newly formed 
stars, the basic feature s seem to be strikingly universal to all de- 
terminations (Kroupa 2001). Initial conditions in star forming 
regions can vary considerably. If the IMF depends on the initial 



conditions, there would thus be no reason for it to be universal. 
Therefore a derivation of the characteristic stellar mass that is 
based on fundamental atomic and molecular physics would be 
more consistent. 

Th ere have been analytical models iTeanslll902l E arson 



19691 IPenstonl Il969t Low & Lvnden-Belll 1 19761 IShdll977t 



Whit worth & Summers 1985) and numerical investigations 
of the effects of various physical processes on collapse 
and fragmentation. These processes include, f or exam- 
ple, magnetic fields JBasu &Mouschoviaslll995t ITomisakal 
Il996t iGalli et alj 1200 lb. feedback from the stars them- 
selves JSilkll995t Nakano et alJl995tfAdams & Fatuzzoll9"96l) 
and competitive coagulation or accretion ( Silk & Takahashi 
19791: iLeieune & Bastienlll986l IPrice & Podsiadlowskil ' 
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Mur ray & Lid 1 19961 iBonnell et alJ 12001 albt burisen et alJ 



2001). In another group of models, initial and environmental 
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conditions, like the structural properties of molecular clouds 



determine the IMF ( Elmegreen & Mathieu 1983; Elmegreen 


1997a b, 


1999, 2000b a, 


2002). 


Larson ( 1973a) and Zinnecker 


( 1984. Il99d 




argued in 


a more statistical approach that the 



central-limit theorem naturally leads to a log-normal stellar 
mass spectrum. Moreover, there are models that connect turbu- 
lent motions in molecular c l ouds to the IMF (e.g. Larson 1981 



| Flecldll982l |Padoan|ll995t IPadoan et allll997t iKlessen et alJ 
ll998[l200rilKlessenl200lUPadoan & Nordlundl2002l) . The in- 
terplay between turbulent motion and self-gravity of the cloud 
leads to a process we call gravoturbulent fragmentation. The 
supersonic turbulence ubiquitously observed in molecular gas 
generates strong density fluctuations with gravity taking over 
in the densest and most massive regions. Once gas clumps be- 
come gravitationally unstable, collapse sets in. The central den- 
sity increases until a protostellar objects form and grow in mass 
via accret ion from the in falli ng envelope. For more detail ed re- 
views see ILarsonl d2003l) and lMac Low & Klessenl d2004l) . 

However, current results are generally based on models that 
do not treat thermal physics in detail. Typically, they use a 
simple equation of state (EOS) which is isothermal with the 
poly tropic exponent y = 1. The true nature of the EOS re- 
mains a major theoretical problem in understanding the frag- 
mentation properties of molecular cl ouds. Some calculations 
invoke cooling during the collapse dMonaghan & Lattanziol 
I199U lTiirneret al.lfl99^ IWhitworth et alJll99.4 . Others in- 
clude radiation transport to account for the heating that oc- 
curs once the cloud reaches den sities of n(Ri) ^ 10 10 cm~ 3 
(Mvhill & Kaul all992llBosll993l) . or simply assume an adia- 
batic equation o f state once that density is ex ceeded (TBonnell 
I1994 iBate etaDll995l) . ISpaans & Silkl d2000t) showed that ra- 
diatively cooling gas can be described by a piecewise poly- 
tropic EOS, in which the polytropic exponent y changes with 
gas density p. Considering a polytropic EOS is still a rather 
crude approximation. In practice the behaviour of y may be 
more complicated and important effects like the temperature of 
the dust, line-trapping and fe edback from newl y-formed stars 
should be taken into account lScaloetal.l ir998). Nevertheless 
a polytropic EOS gives an insight in the differences that a de- 
parture from isotherm ality ev okes. 

Recently iLi et alJ d2003l) conducted a systematic study of 
the effects of a varying polytropic exponent y on gravotur- 
bulent fragmentation. Their results showed that y determines 
how strongly self-gravitating gas fragments. They found that 
the degree of fragmentation decreases with increasing poly- 
tropic exponent y in the range 0.2 < y < 1.4 although the total 
amount of mass in collapsed cores appears to remain roughly 
consistent through this range. These findings suggest that the 
IMF might be quite sensitive to the thermal physics. Earlier, 
one-d imensional simulations bv IPassot & Vazauez-Semadenil 
(1998) already showed that the density probability distribution 
of supersonic turbulent gas displays a dependence on the poly- 
tropic exponent y. However in both computations, y was left 
strictly constant in each case. In this study we extend previous 
work by using a piecewise polytropic equation of state chang- 
ing y at some chosen density. We investigate if a change in y 
determines the characteristic mass of the gas clump spectrum 
and thus, possibly, the turn-over mass of the IMF. 



In Sect.|2]we review what is currently known about the ther- 
mal properties of interstellar gas. In Sect. [3] we approach the 
fragmentation problem analytically, while in Sect. 0] we intro- 
duce our computational method. In Sect. [5] we discuss gravo- 
turbulent fragmentation in non-isothermal gas. In Sect. |S] we 
analyze the resulting mass distribution. We further investigate 
the influence of different turbulent driving fields and different 
scale of driving in Sect.0 Finally, in Sect. [8] we summarize. 

2. Thermal Properties of Star-Forming Clouds 

Gravity in galactic molecular clouds is initially expected to 
be opposed mainly by a combination of superson ic turbu- 
lence and magnetic fields ( Mac Low & Klesserl2004) . The ve- 
locity structure in the clouds i s always observed to be dom- 
inated by large-scale modes (M ac Low & Ossenkopf 2000; 
lOssenkopf et al.ll200lUOssenkopf & Mac Lowll2002i) . In order 
to maintain turbulence for some global dynamical timescales 
and to compensate for gravitational contraction of the cloud as 
a whole, kinetic energy input from external sources seems to be 
required. Star formation then takes place in molecular cloud re- 
gions which are characterized by local dissipation of turbulence 
and loss of magnetic flux, eventually leaving thermal pressure 
as the main force resisting gravity in the small dense prestellar 
cloud cores that actually build u p the stars dKlessen et all2 004: 
IVazquez-Semadeni et alJEoO^) . In agreement with this expec- 
tation, observed prestellar cores typically s how a rough bal- 
ance between gravity an d thermal pressure dBenson & Mversl 
Il989t iMvers etai1ll99lh . Therefore the thermal properties of 
the dense star-forming regions of molecular clouds must play 
an important role in determining how these clouds collapse and 
fragment into stars. 

Early studies of the balance between heating and cooling 
processes in collapsing clouds predicted temperatures of the 
order of 10 K to 20 K, tending to be lower at the higher den- 
sities (e.g.. | Havashi & Nakanolfl965t lHavashilll96rl ILarsonl 
1969, 1973b). In their dynamical collapse calculations, these 
and other authors approximated this somewhat varying tem- 
perature by a simple constant value, usually taken to be 10 K. 
Nearly all subsequent studies of cloud collapse and fragmen- 
tation have used a similar isothermal approximation. However, 
this approximation is actually only a somewhat crude one, valid 
only to a factor of 2, since the temperature is predicted to vary 
by this much above and below the usually assumed constant 
value of 10 K. Given the strong sensitivi ty of t he r esults of 
fragmentation simulations like those of ILi et all d2003l) to the 
assumed equation of state of the gas, temperature variations of 
this magnitude may be important for quantitative predictions of 
stellar masses and the IMF. 

As can be seen in Fig. 2 of ILarsonl d 19851) . observa- 
tional and theoretical studies of the thermal properties of col- 
lapsing clouds both indicate that at densities below about 
10~ 18 gcirr 3 , roughly corresponding to a number density of 
n — 2.5 x 10 s cm~ 3 , the temperature generally decreases with 
increasing density. In this low-density regime, clouds are ex- 
ternally heated by cosmic rays or photoelectric heating, and 
they are cooled mainly by the collisional excitation of low- 
lying levels of C + ions and O atoms; the strong dependence 
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of the cooling rate on density then yields an equilibrium tem- 
pera ture that decrea s es wit h increasing density. The work of 
iKovama & Inutsukal (EoOO). which assumes that photoelectric 
heating dominates, rather than cosmic ray heating as had been 
assumed in earlier work, predicts a very similar trend of de- 
creasing temperature with increasing density at low densi- 
ties. The resulting temperature-density relation can be approxi- 
mated by a power law with an exponent of about -0.275, which 
corresponds to a polytropic equatio n of st ate with y - 0.725. 
The observational results o f lMversN 19781) shown in Fig. 2 of 
ILarsonl i 19851) suggest temperatures rising again toward the 
high end of this low-density regime, but those measurements 
refer mainly to relatively massive and warm cloud cores and not 
to the small, de nse, co l d core s in which low-mass stars form. 
As reviewed by Evans ( 1999), the temperatures of these cores 
are typically only about 8.5 K at a density of 10~ 19 g cm 4 , con- 
sistent with a continuation of the decreasing trend noted above 
and with the continuing validity of a polytropic approximation 
with y « 0.725 up to a density of at least 10~ 19 g cm -3 . 



At densities higher than this, star-forming cloud cores be- 
come opaque to the heating and cooling radiation that de- 
termines their temperatures at lower densities, and at densi- 
ties above 10 _18 gcm~ 3 the gas becomes thermally coupled 
to the dust grains, which then control the temperature by 
their far-infrared thermal emission. In this high-density regime, 
dominated thermally by the dust, there are few direct tem- 
perature measurements because the molecules normally ob- 
served freeze out onto the dust grains, but most of the avail- 
able theoretical predictions are in good agreement c oncern- 
ing the expected thermal behavior of the gas (lLarsonlll9 7 3b ; 
ILow & Lvnden-Belllll976t iMasunaga & Inutsukal 1200(1) . The 
balance between compressional heating and thermal cooling by 
dust results in a temperature that increases slowly with increas- 
ing density, and the resulting temperature-density relation can 
be approximated by a power law with an exponent of about 
0.075, which corresponds to y = 1.075. Taking these values, 
the temperature is predicted to reach a minimum of 5 K at the 
transition between the low-density and the high-density regime 
at about 2 x 10~ 18 gcm ~ 3 , at which p oint the Jeans mass is 
about 0.3 M Q (see also. ILarsonl f2005). The actual minimum 
temperature reached is somewhat uncertain because observa- 
tions have not yet confirmed the predicted very low values, but 
such cold gas would be very difficult to observe; various ef- 
forts to model the observations have suggested central temper- 
atures between 6 K and 10 K for the densest observed prestellar 
cores, whose peak densities may approach 10~ 17 gcm~ 3 (e . g . , 
IZucconiet alJkOOlt lEvans et alJEoOll iTafalla et alJl2004h . A 
power-law approximation to the equation of state with y » 
1.075 is expected to remain valid up to a density of about 
10~ 13 gcrrT 3 , above which increasing opacity to the thermal 
emission from the dust causes the temperature to begin rising 
much more rapidly, resulting in an "opacity limit" on fragmen- 
tation that is somewhat below 0.01 M Q (Low & Lvnden- Bell 
ll976t tMasunag a & Inutsukafe OOO). 



3. Analytical Approach 

Following the above considerations, we use a polytropic equa- 
tion of state to describe the thermal state of the gas in our mod- 
els with a polytropic exponent that changes at a certain critical 
density p c from y\ to y 2 : 



p = k iP * P < Pc 

P = K 2 p 72 p>p c 



(1) 



where K\ and K 2 are constants, and P, and p are thermal pres- 
sure and gas density. For an ideal gas, the equation of state is: 



P = -^- P T(p) 
pm p 



(2) 



where T is the temperature, and k^, P, and m p are Boltzmann 
constant, molecular weight, and proton mass. So the constant 
K can be written as: 



K = -V-r T (p) 
pm v 



(3) 



Since K is defined as a constant in p, it follows for T: 

Ti = flip 71-1 P < Pc 

T 2 = a 2 p^-- 1 p>p c (4) 
where a\ and a 2 are constants. The initial conditions define a\\ 

1-71 



a \ = T()Pq 
At p c it holds that: 

T\(p c ) = r 2 (p c ) 

Thus, a 2 can be written in terms of ci\\ 



(5) 



(6) 



a 2 = aip 7 c ' 72 (7) 

According to the analytical work bv Ijeansi d 19021) on the 
stability of a self-gravitating, isothermal medium the oscilla- 
tion frequency a> and the wavenumber k of small perturbations 
satisfy the dispersion relation 

w 2 - c 2 s k 2 + 4ttGpo = , (8) 

where c s is the sound speed, G the gravitational constant, and po 
the gas density. The perturbation is unstable if the wavelength A 
exceeds the Jeans length /lj = 2n/kj or, equivalently, if the mass 
exceeds the Jeans mass 



An I A 

Ms = yp 



, ) ,3 _5/2 



(9) 



Note that we define the Jeans mass Mj as the mass originally 
contained within a sphere of diameter As. 

In a system with a polytropic EOS, i.e., P = Kp 7 , the sound 
speed is 



dP 
dp 



1/2 



1/2 r r -i)/2 



Thus, the Jeans mass can be written as 

n V2 /K ,V2 



Ms 



7!""" IKV 



y 3/2 p (3/2) y -2 



(10) 



(11) 
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Fig. 1. Local Jeans mass as a function of density for four runs 
with different critical densities « c . For comparison the depen- 
dence is also shown for the isothermal case {dotted line). The 
Jeans mass changes at the critical density. The initial mean den- 
sity and the density at which sink particles form are represented 
by the vertical dashed lines. The dashed-dotted lines show the 
minimal resolvable mass for the runs with the highest resolu- 
tion. 



Using Eqs. 3, 4, 5, 7 and 1 1 one finds: 



My 



M J2 = 




-,3/2(3/2)7, -2 



p <p c 



r 3/2 p (3/2)( yi - y2 ) p( 3/ 2)y2 -2 p>pc 



The sound speed changes when the polytropic index changes at 
p c , so Mj also varies (see Fig.Q, such that: 



M J2 



V72 



(12) 



If we use 7i=0.7 and 72-1-1 as justified in Sect.[5]then M n oc 

p-°- 95 . 

During the initial phase of collapse, the turbulent flow pro- 
duces strong ram pressure gradients that form density enhance- 
ments. Higher density leads to smaller local Jeans masses, so 
these regions begin to collapse and fragment. Simulations with 
an SPH code different from the one used in the present work 
show that fragmentation occurs more efficiently for smaller 
values of y, and less effici ently for y > 1, cutting o ff en- 
tirely at 7 > 1.4 dli et alJ EooH lArcoragi et alJ fl99 1). For 
filamentary systems, fragm entation already stops for y > 1 
jKawachi & Hanawal 19981) . This point is discussed in more de- 
tail in Sec.lU 



Table 1. Sample parameters, name of the environment used in 
the text, driving scale k, critical density n c , number of SPH par- 
ticles, number N of protostellar objects (i.e., "sink particles" in 
the centers of protostellar cores) at final stage of the simulation, 
percentage of accreted mass at final stage M acc /M tot 



Nfttnc 


k 


cm 


nartinlp 

Lfui t 1 ^ 1 V- 

number 


N 


M. icc 
*>lol 

|%J 


Ik2 


1..2 




205 379 


59 


56 


Ik2b 


1..2 




1 000 000 


73 


78 


Ik2L 


1..2 




9 938 375 


6 


4 


R5k2 


1..2 


4.63 


205 379 


22 


73 


R5k2b 


1..2 


4.63 


1 000 000 


22 


70 


R6k2 


1..2 


5.63 


205 379 


64 


93 


R6k2b 


1..2 


5.63 


1 000 000 


54 


61 


R7k2 


1..2 


6.63 


205 379 


122 


84 


R7k2b 


1..2 


6.63 


1 000 000 


131 


72 


R7k2L 


1..2 


6.63 


1 953 125 


143 


46 


R8k2 


1..2 


7.63 


205 379 


194 


78 


R8k2b 


1..2 


7.63 


1 000 000 


234 


53 


R8k2L 


1..2 


7.63 


5 177 717 


309 


29 


R5k8 


7.. 8 


4.63 


205 379 


1 


64 


R6k8 


7.. 8 


5.63 


205 379 


38 


68 


R7k8 


7. .8 


6.63 


205 379 


99 


60 


R8k8 


7.. 8 


7.63 


205 379 


118 


72 


R5k2rl 


1..2 


4.63 


205 379 


16 


62 


R6k2rl 


1..2 


5.63 


205 379 


34 


72 


R7k2rl 


1..2 


6.63 


205 379 


111 


68 


R8k2rl 


1..2 


7.63 


205 379 


149 


64 


R5k2r2 


1..2 


4.63 


205 379 


21 


72 


R6k2r2 


1..2 


5.63 


205 379 


51 


74 


R7k2r2 


1..2 


6.63 


205 379 


119 


70 


R8k2r2 


1..2 


7.63 


205 379 


184 


70 


R5k2r3 


1..2 


4.63 


205 379 


18 


90 


R6k2r3 


1..2 


5.63 


205 379 


52 


85 


R7k2r3 


1..2 


6.63 


205 379 


123 


76 


R8k2r3 


1..2 


7.63 


205 379 


196 


71 



What happens when y increases above unity at the critical 
density p c ? One suggestion is that the increase in y is sufficient 
to strongly reduce fragmentation at higher densities, introduc- 
ing a characteristic scale into the mass spectrum at the value 
of the Jeans mass at p c . Then the behavior of the Jeans mass 
with increasing critical density would immediately allow us to 
derive the scaling law 

-0.95 



M ch oc p c 



(13) 



This simple analytical consideration would then predict a 
characteristic mass scale which corresponds to a peak of the 
IMF at 0.35 M for a critical density of p c = 10~ 18 g cirT 3 or 
equivalently a number density of « c = 2.5 x 10 5 cirT 3 when 
using a mean molecular weight yi = 2.36 appropriate for solar 
metallicity molecular clouds in the Milky Way. Note, however, 
that this simple scaling law does not take any further dynamical 
processes into account. 

4. Numerical Method 

In order to test the scaling of the characteristic mass M c h given 
by Eq. [O] we carry out simulations of regions in turbulent 
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Fig. 2. Temperature as a function of density for four runs with 
different critical densities n c . The dotted lines show the ini- 
tial conditions. The curve has a discontinuous derivative at the 
critical density. 



molecular clouds in which we vary the critical density p c and 
determine the resulting mass spectra of protostellar objects, and 
thus M<±. During gravoturbulent fragmentation it is necessary 
to follow the gas over several orders of magnitude in density. 
The method of choice therefore is smoothed particle hydrody- 
namics (SPH). Excellent overviews of the method, its numer- 
ical imple mentation, an d s ome of its a pplications are given in 
reviews bv lBend Jl990l) and lMona ghanM 1 9921). We use the par- 
allel code GADGET, designed bv lSpringel et aDfeOOll) . SPH is 
a Lagrangian method, where the fluid is represented by an en- 
semble of particles, and flow quantities are obtained by averag- 
ing over an appropriate subset of SPH particles. We use a spher- 
ically symme tric cubic splin e function to define the smooth- 
ing kernel (Monag hanlll992l) . The smoothing length can vary 
in space and time, such that the number of considered neigh- 
bors is always approximately 40. The method is able to resolve 
large density contrasts as particles are free to move, and so 
naturally the particl e concentration i ncreas es in high-density 
regions. We use the Bate & Burkert ( 1997) criterion to deter- 
mine the resolution limit of our calculations. It is adequate for 
the problem considered here, where we follow the evolution of 
highly nonlinear density fluctuations created by supersonic tur- 
bulence. We have performed a resolution study with up to 10 7 
SPH particles to confirm this result. 



4.1. Sink particles 

SPH simulations of collapsing regions become slower as more 
particles move to higher density regions and hence have small 
timesteps. Replacing dense cores by one single particle leads 
to considerable increase of the overall computational perfor- 



mance. Introducing sink particles allows us to follow the dy- 
namical evolution of the system over many free-fall times. 
Once the density contrast in the center of a collapsing cloud 
core exceeds a value of 5000, corresponding to a threshold den- 
sity of about % = 4x 10 8 cm" 3 (see Fig.[T), th e entire ce ntral 
region is replaced by a "sink particle" (B ate. Bonnell. & Pried 
Il995h . It is a single, non-gaseous, massive particle that only in- 
teracts with normal SPH particles via gravity. Gas particles that 
come within a certain radius of the sink particle, the accretion 
radius r acc , are accreted if they are bound to the sink particle. 
This allows us to keep track of the total mass, the linear and 
angular momentum of the collapsing gas. 

Each sink particle defines a control volume with a fixed 
radius of 310 AU. This radius is chosen such that we always 
resolve the Jeans scale belo w the threshold density n t h, fol- 
lowing [Bat^^iuricertl Jl997l) . We cannot resolve the subse- 
quent evolution in its interior. Combination with a detailed one- 
dimensional implicit radiation hydrodynamic method shows 
that a protosta r forms in the very ce nter about 10 3 yr after sink 
creation ( Wuchterl & Klessen 2001). We subsequently call the 
sink protostellar object or simply protostar. Altogether, the 
sink particle represents only the innermost, highest-density part 
of a larger collapsing region. The technical details on the imple- 
mentation of sink particles in the parallel SPH code GADGET 
can be found in Appendix A. 

Protostellar collapse is accompanied by a substantial loss 
of spe cific angular momentum, e ven in the absence of magnetic 
fields (Jappsen & Klessen 12004 . Still, most of the matter that 
falls in will assemble in a protostellar disk. It is then transported 
inward by viscous and possibly gravitational torques (e.g.. 



Bodenheimer 1995; Papaloizou & Lin 1995: lLin&Papaloizoul 
1996). With typical disk sizes of order of several hundred AU, 
the control volume fully encloses both star and disk. If low an- 
gular momentum material is accreted, the disk is stable and 
most of the material ends up in the central star. In this case, 
the disk simply acts as a buffer and smooths eventual accre- 
tion spikes. It will not delay or prevent the mass growth of 
the central star by much. However, if material that falls into 
the control volume carries large specific angular momentum, 
then the mass load onto the disk is likely to be faster than in- 
ward transport. The disk grows large and may become gravita- 
tionally unstable and fragment. This m ay lead to the f ormation 
of a binary or higher -order multiple (Bod enheimer et all f2000: 
iFromang et all2002l) . 

4.2. Model parameters 

We include turbulence in our version of t he code that is 
driven unif ormly with t he me thod described bv lMac Low et alj 
Jl998t) and iMac Low! Jl999l) . The observed turbulent veloc- 
ity field in molecular clouds wil l decay in a crossing time 
if not continuously replenished (|M ac Low & Klessenf l2004 ; 
IScalo & Elmegreenl2004HElmegreen & Scalol2004 . If the tur- 
bulence decays, then only thermal pressure prevents global col- 
lapse. In our case, we examine regions globally supported by 
the turbulence at the initial time. We choose to model con- 
tinuously driven turbulence leading to inefficient star forma- 
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Fig. 3. Column density distribution of the gas and location of identified protostellar objects (black circles) using the high- 
resolution models R6..8k2b at the stage where approximately 50% of the gas is accreted. Projections in the xy-, xz-, and yz-plane 
are shown for three different critical densities. 



tion rather than a globally collaps ing region producing efficient 
star formation jMac Low & Klessenll2004l) . This is achieved 
here by applying a nonlocal driving schem e that inserts en- 
ergy in a limited range of wave numbers fe. lMac Low! ill 9991) 
shows that hydrodynamical turbulence decays with a constant 
energy-dissipation coefficient. Thus a constant kinetic energy 
input rate E m is able to maintain the observed turbulence and to 
stabilize the system against gravitational contraction on global 
scales. The driving strength is adjusted to yield a constant tur- 
bulent rms Mach number Alrms = 3.2. Furthermore, it is known 
that the velocity structure in molecula r clouds is always domi- 
nated by the largest-scale mo d e (e.g-.lMac Low & Ossenkopfl 
2000 t lOssenkopfetalJ 1200 it lOssenkopf & Mac Low! I2002T 
Brunt et al 1120041) . consequently we insert energy on scales of 



the order of the size of our computational domain, i.e. with 
wave numbers k = 1..2. In the adopted integration scheme, we 
add the turbulent energy AEm = E[ n At to the system at each 
timestep At whenever more than 40% of the SPH particles are 
advanced. 



In all our models we adopt an initial temperature of 1 1 .4 K 
corresponding to a sound speed c s = 0.2kms _1 , a molecu- 
lar weight /j of 2.36 and an initial number density of n — 
8.4 x 10 4 crrr 3 , which is typical for star-forming molecular 
cloud regions (e.g., p-Ophiuchi, see Motte et al. 1998 or the 
central region of the Orion Nebula Cluster, see Hillenbrand 
1997; Hillenbrand & Hartmann 1998). Our simulation cube 
holds a mass of 120 M Q and has a size of L = 0.29 pc. The 



A.-K. Jappsen et al.: Mass spectrum from non-isothermal gravoturbulent fragmentation 



7 



cube is subject to periodic boundary conditions in every direc- 
tion. The mean initial Jeans mass is (Mj); = 0.7 M . 

We use the EOS described in Sect. [5] and compute models 
with 4.3 Xl0 4 cm 3 < n c < 4.3 xlO 7 cm 3 . Note, that the lowest 
and the highest of these critical densities represent rather ex- 
treme cases. From Fig. [2] where we show the temperature as a 
function of number density, it is evident that they result in tem- 
peratures that are too high or too low compared to observations 
and theoretical predictions. Nevertheless, including these cases 
clarifies the existence of a real trend in the dependence of the 
characteristic mass scale on the critical density. Each simula- 
tion starts with a uniform density. Driving begins immediately, 
while self-gravity is turned on at t = 2.0 fff , after turbulence is 
fully established. The global free-fall timescale is tg » 10 5 yr. 
Our models are named mnemonically. R5 up to R8 stand for 
the critical density n c (4.3 x 10 4 cirT 3 < n c < 4.3 x 10 7 cm~ 3 ) 
in the equation of state, k2 or k8 stand for the wave numbers 
(k = 1..2 or k — 7. .8) at which the driving energy is injected 
into the system and b flags the runs with 1 million gas particles. 
The letter L marks the high resolution runs for critical densities 
h c = 4.3 x 10 6 cm 4 and n c = 4.3 x 10 7 cm 4 with 2 million and 
5.2 million particles, respectively. Different realizations of the 
turbulent velocity field are denoted by rl, r2, r3. For compari- 
son we also run isothermal simulations marked with the letter I 
that have particle numbers of approximately 200000, 1 million 
and 10 million gas particles. 

The number of particles determines the minimal resolv- 
able Jeans mass in our models. Figure^shows the dependence 
of the local Jeans mass on the density for the adopted poly- 
tropic equation of state. At the critical density the dependency 
of the Jeans mass on density changes its behavior. The mini- 
mum Jeans mass M res that needs to be resolved occurs at the 
density at which sink particles are formed. A local Jeans mass 
is considered resolved if it contains at least 2 x A^ ne i g h = 80 
SPH particles (Bat e & Burkerll997h . As can be seen in Fig.[T] 
we are able to resolve M les with 1 million particles for critical 
densities up to n c = 4.3 x 10 5 cirT 3 . Since this is not the case 
for n c = 4.3 x 10 6 cirT 3 and n c = 4.3 x 10 7 cirT 3 , we repeat 
these simulations with 2 million and 5.2 million particles, re- 
spectively. Due to long calculation times we follow the latter 
only to the point in time where about 30% of the gas has been 
accreted. 

At the density where y changes from below unity to above 
unity, the temperature reaches a minimum. This is reflected 
in the "V"-shape shown in Figure [2] All our simulations start 
with the same initial conditions in temperature and density as 
marked by the dotted lines. 

In a further set of simulations we analyze the influence of 
changing the turbulent driving scheme on fragmentation while 
using a polytropic equation of state. These models contain 2 x 
10 5 particles e ach. 

Following Bate & Burkert ( 1997), these runs are not con- 
sidered fully resolved for « c > 4.3 x 10 5 at the density of 
sink particle creation, since M les falls below the critical mass 
of 80 SPH particles. We note, however, that the global accre- 
tion history is not strongly affected. First, we study the effect 
of different realizations of the turbulent driving fields on typ- 
ical masses of protostellar objects. We simply select different 
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Fig. 4. Temporal evolution of the number of protostellar ob- 
jects (upper plot ) and of the ratio of accreted gas mass to total 
gas mass (lower plot ) for models R5..8k2b. The legend shows 
the logarithms of the respective number densities in crrT 3 . The 
times are given in units of a free-fall time Ts. We also show 
the models R7k2L and R8k2L with 2 x 10 6 and 5 x 10 6 par- 
ticles, respectively, which are denoted by the letter "L". For 
comparison the dotted lines indicate the values for the isother- 
mal model Ik2b. 

random numbers to generate the field while keeping the over- 
all statistical properties the same. This allows us to assess the 
statistical reliability of our results. These models are labeled 
from R5..8k2rl to R5..8k2r3. Second, driving in two different 
wavenumber ranges is considered. Most models are driven on 
large scales (1 < k < 2) but we have run a set of models driven 
on small scales (7 < k < 8) for comparison. 

The main model parameters are summarized in Tabled 

5. Gravoturbulent fragmentation in polytropic gas 

Turbulence establishes a complex network of interacting 
shocks, where converging flows and shear generate filaments of 
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high density. The interplay between gravity and thermal pres- 
sure determines the further dynamics of the gas. Adopting a 
polytropic EOS (Eq. [Q, the choice of the polytropic exponent 
plays an important role determining the fragmentation behav- 
ior. From Eq.^Jit is evident that y = 4/3 constitutes a critical 
value. A Jeans mass analysis shows that for three-dimensional 
structures Mj increases with increasing density if y is above 
4/3. Thus, y > 4/3 results in the termination of any gravita- 
tional collapse. 

Also, collapse and fragmentation in filaments depend on 
the equation of state. The equilibrium and stability of fil- 
ame ntary structures has been studie d extensively, beginning 
with Chandrasekhar & Fermil dl953h . and this work has been 
reviewed by iLarsonl Q1985 , 2003J). For many types of col- 
lapse problems, insight into the qualitative behavior of a col- 
lapsing config uration can be gained from similarity solutions 
(Larson 2003). For the collapse of cylinders with an assumed 
polytropic equation of stat e solutions have been derived by 
iKawachi & Hanawal dl 998h . and these authors found that the 
existence of such solutions depends on the assumed value of y: 
similarity solutions exist for y < 1 but not for y > 1. These 
authors also found that for y < 1, the collapse becomes slower 
and slower as y approaches unity from below, asymptotically 
coming to a halt when y — 1 . This result shows in a particu- 
larly clear way that y — 1 is a critica l case for the collapse of 
filaments. Kawachi & Hanawa ( 1998) suggested that the slow 
collapse that is predicted to occur for y approaching unity will 
in reality cause a filament to fragment into clumps, because 
the timescale for fragmentation then becomes shorter than the 
timescale for collapse toward the axis of an ideal filament. If the 
effective value of y increases with increasing density as the col- 
lapse proceeds, as is expected from the predicted thermal be- 
havior discussed in Sect.|2] fragmentation may then be particu- 
larly favored to occur at the densi ty whe re y approaches unity. 
In their numerical study ILi et al.l ll2003h found, for a range of 
assumed polytropic equations of state, that the amount of frag- 
mentation that occurs is indeed very sensitive to the value of 
the polyt ropic exponent j, esp ecially for values of y near unity 
(see also. lArcoragi et alJl99ll) . 

The fact that filamentary structure is so prominent in our 
results and other simulations of star formation, together with 
the fact that most of the stars in these simulations form in fil- 
aments, suggests that the formation and fragmentation of fila- 
ments may be an important mode of star formation quite gen- 
erally. This is also supported by the fact that many observed 
star-forming clouds have filamentary structure, and by the evi- 
dence that much of the star formation in these clouds occurs in 
filaments ( Schneider & Elmeareen 1979; Larson 1985; Curry 
l2002llHartmanr|2002h As w e note in Sects. |2]and|3]and fol- 
lowing lLarsonl ill973hi r2005). the Jeans mass at the density 
where the temperature reaches a minimum (see Fig. |2}, and 
hence, y approaches unity, is predicted to be about 0.3 M Q , co- 
incidentally close to the mass at which the stellar IMF peaks. 
This similarity is an additional hint that filament fragmentation 
with a varying polytropic exponent may play an important role 
in the origin of the stellar IMF and the characteristic stellar 
mass. 



The filamentary structure that occurs in our simulations is 
visualized in Fig. [3] Here we show the column density distribu- 
tion of the gas and the distribution of protostellar objects. We 
display the results for three different critical densities in xy-, 
xz- and yz-projection. The volume density is computed from 
the SPH kernel in 3D and then projected along the three prin- 
cipal axes. Figure |3 shows for all three cases a remarkably fil- 
amentary structure. These filaments define the loci where most 
protostellar objects form. 

Clearly, the change of the polytropic exponent y at a certain 
critical density influences the number of protostellar objects. 
If the critical density increases then more protostellar objects 
form but the mean mass decreases. 

We show this quantitatively in Fig. |4] In (a) we compare 
the number of protostellar objects for different critical densi- 
ties n c for the models R5...8k2b. The rate at which new proto- 
stars form changes with different n c . Models that switch from 
low y to high y at low densities built up protostellar objects 
more rarely than models that change y at higher densities. 

Figure|4]3 shows the accretion histories (the time evolution 
of the combined mass fraction of all protostellar objects) for 
the models R5..8k2b. Accretion starts for all but one case ap- 
proximately at the same time. In model R5k2b, y - 1.1 already 
at the mean initial density, thus y does not change during col- 
lapse. In this case accretion s tarts at a later time. This confirms 
the finding bv lLi et all £2003 ) that accretion is delayed for large 
y. In the other four cases the accretion history is very similar 
and the slope is approximately the same for all models. 

In both plots we also show the results from our high resolu- 
tion runs R7k2L and R8k2L. These simulations with 2 million 
and 5.2 million particles, respectively, have an accretion his- 
tory similar to the time evolution of the accreted mass fraction 
in the runs with 1 million particles. The number of protostellar 
objects, however, is larger for the runs with increased particle 
numbers. Combining our results in these two figures we find 
that an environment where y changes at higher densities pro- 
duces more, but less massive objects. Thus, the mean mass of 
protostellar objects does indeed depend on the critical density 
where y changes from 0.7 to 1.1. 

6. Dependence of the characteristic mass on the 
equation of state 

Further insight into how the characteristic stellar mass may 
depend on the critical density can be gained from the mass 
spectra of the protostellar objects, which we show in Fig. [5] 
We plot the mass spectra of models R5...6k2b, model R7k2L 
and model R8k2L at different times, when the fraction of mass 
accumulated in protostellar objects has reached approximately 
10%, 30% and 50%. In the top row we also display the results 
of an isothermal run for comparison. We used the same initial 
conditions and parameters in all models shown. Dashed lines 
indicate the specific mass resolution limits. 

W e find closest correspondence with the ob served IMF 
(see, IScalolll998t iKroupalEool IChabrieJEool for a crit- 
ical density of 4.3 x 10 6 cm~ 3 and for stages of accretion 
around 30% and above. At high masses, our distribution fol- 
lows a Salpeter-like power law. For comparison we indicate 
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Fig. 5. Mass spectra of protostellar objects for models R5..6k2b, model R7k2L and model R8k2L at 10%, 30% and 50% of total 
mass accreted on these protostellar objects. For comparison we also show in the first row the mass spectra of the isothermal 
run Ik2b. Critical density n c , ratio of accreted gas mass to total gas mass M acc /M tot and number of protostellar objects are given 
in the plots. The vertical solid line shows th e position of the median mass. The dotted line has a slope of -1.3 and serves as a 
reference to the Salpeter value (Salp eteJl95i . The dashed line indicates the mass resolution limit. 



the Salpeter slope x « 1. 3 iSalr)eterlll95.i where the IMF is 
defined by dA^/dlogm oc m~ x . For masses about the median 
mass the distribution exhibits a small plateau and then falls off 
towards smaller masses. 



The model R5k2b where the change in y occurs below 
the initial mean density, shows a flat distribution with only 
few, but massive protostellar objects. They reach masses up 
to 10 M and the minimal mass is about 0.3 M Q . All other 
models build up a power-law tail towards high masses. This is 
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due to protostellar accretion processe s, as more and more gas 
gets turned into stars (see also. iBonnell et al1l2001bt iKlessenl 
1200 ll LSchmeia & Klessenl 120041) ." The distribution becomes 
more peaked for higher n c and there is a shift to lower masses. 
This is already visible in the mass spectra when the protostellar 
objects have only accreted 10% of the total mass. Model R8k2L 
has minimal and maximal masses of 0.013 M Q and 1.0 M , re- 
spectively. There is a gradual shift in the median mass (as in- 
dicated by the vertical line) from Model R5k2b, with M me d = 
2.5 M at 50%, to Model R8k2L, with M med = 0.05 M Q at 50%. 
A similar trend is noticeable during all phases of the model evo- 
lution. 

This change of median mass with critical density « c is 
depicted in Fig. |6] Again, we consider models R5...6k2b, 
model R7k2L and model R8k2L. The median mass decreases 
clearly with increasing critical density as expected. As we re- 
solve higher density contrasts the median collapse mass de- 
creases. We fit our data with a straight line. The slope takes 
values between -0.4 and -0.6. These values are larger than the 
slope -0.95 derived from the simple scaling law (Eq.ll3> based 
on calculating the Jeans mass Mj at the critical density « c 

One possible reason for this deviation is the fact that most 
of the protostellar objects are members of tight groups. Hence, 
they are subject to mutual interactions and competitive accre- 
tion that may change the environmental context for individual 
protostars. This in turn i nfluences the mass distr ibution and its 
characteristics (see, e.g., IBonnell et al.l2001albl) . Another pos- 
sible reason is that the mass that goes into filaments and then 
into collapse may depend on further environmental parameters, 
some of which we discuss in Sect. [7] 

7. Dependence of the characteristic mass on 
environmental parameters 

7.1. Dependence on realization of the turbulent 
velocity field 

We compare models with different realizations of the turbulent 
driving field in Fig. We fit our data with straight lines for 
each stage of accretion. Figure shows the results of mod- 
els R5..8k2 which were calculated with the same parameters 
but lower resolution than the models used for Fig. [6] As dis- 
cussed in Sec. [5] although the number of protostellar objects 
changes with the number of particles in the simulation, the 
time evolution of the total mass accreted on all protostellar 
objects remains similar. Thus, lower resolution models exhibit 
the same general trend as their high-resolution counterparts and 
show the same global dependencies. We notice, however, that 
the slope of the M me dian-»c relation typically is shallower in 
the low-resolution models. This can be seen when comparing 
Fig- Eland Fig.f^, where we used identical turbulent driving 
fields. 

For all four different realizations of the turbulent driving 
field shown in Fig.f^-d we see a clear trend of decreasing me- 
dian mass M me dian with increasing n c . We conclude that, qual- 
itatively, the M me dian-«c relation is independent of the details 
of the turbulent driving field but, quantitatively, there are sig- 
nificant variations. This is not surprising given the stochastic 
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Fig. 6. The plot shows the median mass of the protostellar ob- 
jects over critical density for models R5..6k2b, model R7k2L 
and model R8k2L. We display results for different ratios of ac- 
creted gas mass to total gas mass M acc /M tot . We fit our data with 
straight lines for different stages of accretion. The slopes have 
the following values: -0.43 + 0.05 (solid line), -0.52 + 0.06 
(dashed-dotted line), -0.60 + 0.07 (dashed line). 

nature of tu rbulent flows. A furth er d iscussion on this issu e can 
be found in lKlessen etall kOOdti and lHeitsch etall fcOOll) . 

7.2. Dependence on the scale of turbulent driving 

Figure^ shows the results for models where the driving scale 
has been changed to a lower value (k = 7..8). The overall 
dependency of M me dian on n c is very similar to the cases of 
large-scale turbulence. However, we note considerably larger 
uncertainties in the exact value of the slope. This holds es- 
pecially for the phases where 30% and 50% of the total gas 
mass has been converted into stars. One of the reasons is lower 
statistics, i.e., for the smallest critical density only one pro- 
toste llar object fo rms. Moreover, it has already been noted 
by iLi et al.l (12003) that driving on small wavelength results 
in less fragmentation. The forming small-scale density struc- 
ture is not so strongly filamentary, compared to the case of 
small-scale driven turbulence. Local differences have a larger 
influence on the results for driving on small wavelengths. 
Nevertheless, for most of the models the mean mass decreases 
with increasing critical density. Observational evidence sug- 
gest s that real molecular clouds are driven from la rge scales 
(e.g.lOssenkonf & Mac LowfeoTiilBmnt et akEool. 

8. Summary 

Using SPH simulations we investigate the influence of a piece- 
wise poly tropic EOS on fragmentation of molecular clouds. We 
study the case where the polytropic index y changes from a 
value below unity to one above at a critical density n c . We con- 
sider a broad range of values of n c around the realistic value to 
determine the dependence of the mass spectrum on n c . 

Observational evidence predicts that dense prestellar cloud 
cores show rough balance between gravity and thermal pres- 
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sure dBenson & Mversl Il989t iMvers et alJ Il99ll) . Thus, the 
thermodynamical properties of the gas play an important role 
in determining how dense star-forming regions in molecular 
clouds collapse and fragment. Observational and theoretical 
studies of the thermal properties of collapsing clouds both 
indicate that at densities below lCT 18 gcirT 3 , roughly corre- 
sponding to a number density of n c = 2.5 x 10 5 cnT 3 , the 
temperature decreases with increasing density. This is due to 
the strong dependence of m olecular cooling rates on density 
iKovama & Inutsukal EoOO). Therefore, the polytropic expo- 
nent y is below unity in this density regime. At densities 
above 10~ 18 gcm~ 3 , the gas becomes thermally coupled to 
the dust grains, which then control the temperature by far- 
infrared thermal emission. The balance between compressional 
heating and thermal cooling by dust causes the temperature 
to increase again slowly with increasing density. Thus the 
tempe rature-d ensity relation can be approximated with y above 
unity jLarsonlll985l) in this regime. Changing y from a value 
below unity to a value above unity results in a minimum tem- 
perature at the critical density. ILi et al.l J2003ll showed that gas 
fragments efficiently for y < 1 .0 and less efficiently for higher 
y. Thus, the Jeans mass at the critical density defines a char- 
acteristic mass for fragmentation, which may be related to the 
peak of the IMF. 

We investigate this relation numerically by changing y 
from 0.7 to 1.1 at different critical densities n c varying from 
4.3 x 10 4 cm~ 3 to 4.3 x 10 7 cm~ 3 . A simple scaling argument 
based on the Jeans mass Mj at the critical density n c leads to 
M c h oc n~ . If there is a close relation between the average 
Jeans mass and the characteristic mass of a fragment, a similar 
relation should hold for the expected peak of the mass spec- 
trum. Our simulations qualitatively support this hypothesis, 
however, with the weaker density dependency M c h oc n~° ■ 5±0 - 1 , 
The density at which y changes from below unity to above 
unity selects a characteristic mass scale. Consequently, the 
peak of the resulting mass spectrum decreases with increasing 
critical density. This spectrum not only shows a pronounced 
peak but also a powerlaw tail towards higher masses. Its be- 
havior is thus similar to the observed IMF. 

Altogether, supersonic turbulence in self-gravitating 
molecular gas generates a complex network of interacting 
filaments. The overall density distribution is highly inho- 
mogeneous. Turbulent compression sweeps up gas in some 
parts of the cloud, but other regions become rarefied. The 
fragmentation behavior of the cloud and its ability to form stars 
depend strongly on the EOS. If collapse sets in, the final mass 
of a fragment depends not only on the local Jeans criterion, but 
also on additional processes. For example, protostars grow in 
mass by accretion from their surrounding material. In turbulent 
clouds the properties of the gas reservoir are continuously 
changing. In a dense cluster environment, furthermore, proto- 
stars may interact with each other, leading to ejection or mass 
exchange. These dynamical factors modify the resulting mass 
spectrum, and may explain why the characteristic stellar mass 
depends on the EOS more weakly than expected. 

We also studied the effects of different turbulent driving 
fields and of a smaller driving scale. For different realizations 
of statistically identical large-scale turbulent velocity fields we 



consistently find the characteristic mass is decreasing with in- 
creasing critical mass. However, there are considerable varia- 
tions. The influence of the natural stochastic fluctuations in the 
turbulent flow on the resulting median mass is almost as pro- 
nounced as the changes of the thermal properties of the gas. 
Also when inserting turbulent energy at small wavelengths we 
see the peak of the mass spectrum decreasing with increasing 
critical density. 

Our investigation supports the idea that the distribution of 
stellar masses depends, at least in part, on the thermodynamic 
state of the star-forming gas. If there is a low-density regime 
in molecular clouds where temperature T sinks with increas- 
ing density p, followed by a higher-density phase where T in- 
creases with p, fragmentation seems likely to be favored at the 
transition density where the temperature reaches a minimum. 
This defines a characteristic mass scale. The thermodynamic 
state of interstellar gas is a result of the balance between heat- 
ing and cooling processes, which in turn are determined by fun- 
damental atomic and molecular physics and by chemical abun- 
dances. The derivation of a characteristic stellar mass can thus 
be based on quantities and constants that depend solely on the 
chemical abundances in a molecular cloud. The current study 
using a piecewise polytropic EOS can only serve as a first step. 
Future work will need to consider a realistic chemical network 
and radiation transfer processes in gas of varying abundances. 
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Appendix A: Implementation of Sink Particles 

The parallel version of GADGET distributes the SPH particles onto 
the individual processors, using a spatial domain decomposition. 
Thus, each processor hosts a rectangular piece of computational vol- 
ume. If the position of a sink particle is near the boundary of this 
volume, the accretion radius overlaps with domains on other proces- 
sors. We therefore communicate the data of the sink to all processors. 
Every processor searches for gas particles within the accretion radius 
of the sink. Three criteria determine whether the particle gets accreted 
or not. First, the particle must be bound to the sink particle, i.e., the 
kinetic energy must be less than the magnitude of the gravitational 
energy. Second, the specific angular momentum of the particle must 
be less than what is required to move on a circular orbit with radius 
r mc around the sink particle. Finally, the particle must be more tightly 
bound to the candidate sink particle than to other sink particles. Once 
the central region of a collapsing gas clump exceeds a density contrast 
ZS.pl p ~ 5000, we introduce a new sink particle. The procedure for dy- 
namically creating a sink particle is as follows. We search all proces- 
sors for the gas particle with the highest density. When this density is 
above the threshold and when its smoothing length is less than half the 
accretion radius, then the gas particle is considered to become a sink 
particle. If the accretion radius around the candidate particle overlaps 
with another domain, its position is sent to the other processors. Every 
processor searches for the particles that exist in its domain and, simul- 
taneously, within the accretion radius of the candidate particle. These 
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particles and the candidate particle undergo a series of tests to decide 
if they should form a sink particle. First, the new sink particle must be 
the only one within two accretion radii. Second, the ratio of thermal 
energy to the magnitude of the gravitational energy must be less than 
0.5. Third, we require that the total energy is less than zero. Finally, 
the divergence of the accelerations on the particles must be less than 
zero. If all these tests are passed, the particle with the highest density 
turns into a sink particle with position, velocity and acceleration de- 
rived from the center of mass values of the original gas particles. If 
these original particles are distributed over several processors the cen- 
ter of mass values have to be communicated correctly to the processor 
that hosts the new sink particle. 

Ideally, the creation of sink particles in an SPH simulation should 
not affect the evolution of the gas outside its accretion radius. In prac- 
tice there is the discontinuity in the SPH particle distribution due to 
the hole produced by the sink particles. This affects the pressure and 
viscous forces on particles outside. We have implemented adequate 
boundary conditions at the "surface" of the sink particles as described 
in detail in lBateetalJil995T) to correct for these effects. 

Fo llowing bate et all ll995t) we use the iBoss & Bodenheimed 
ll979l) standard isothermal test case for the collapse and fragmentation 
of an interstellar cloud core to check our implementation. Initially, the 
cloud core is spherically symmetric with a small m = 2 perturbation 
and uniformly rotating. As gravitational collapse proceeds a rotation- 
ally supported high-density bar builds up in the center embedded in 
a disk-like structure. The two ends of the bar become gravitationally 
unstable, resulting in the formation of a binar y system . We see no 
further subfragmentation (see also. iTruelove et al]|l997l) . These tests 
show that the precise creation time and the mass of the sink particle 
at the time of its formation can vary somewhat with the number of 
used processors. We also find that simulations with different proces- 
sor numbers show small deviations in the exact positions and veloc- 
ities of the gas particles. These variations are due to the differences 
in the extent of the domain on each processor. When the force on a 
particular particle is computed, the force exerted by distant groups 
is approximated by their lowest multipole moments. Since each pro- 
cessor constructs its own Barnes and Hut tree differences in the tree 
walk result in differences in the computed force. Hence, the formation 
mass and time of sink particles depend on the computational setup. 
Nevertheless, these differences are only at the 0.1% level and the to- 
tal number of collapsing objects is not influenced by a change in the 
number of processors. 
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Fig. 7. Median mass of protostellar objects over critical density 
at different evolutionary phases (ratio of accreted gas mass over 
total gas mass M acc /M tot :10% (circle), 30% (star), 50% (trian- 
gle), 70% (square)). In (a)-(d) we show identical models but 
with different turbulent velocity fields (same root-mean-square 
velocity). The model in (e) is driven on a smaller scale k = 7. .8 
than in the other cases. All relevant parameters are summarized 
in Tabled We fit the median values with straight lines for dif- 
ferent stages of accretion. The slopes are given in the plot and 
denoted as follows: 10% - solid line, 30% - dashed-dotted line, 
50% - dashed line, 70% - dashed-double-dotted line. 



